home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / CUGUK / C005.ZIP / ATAN.C < prev    next >
Text File  |  1990-01-19  |  3KB  |  117 lines
  1. /********************************************************************
  2.  * C Users Group (U.K) C Source Code Library File CUGLIB.005        *
  3.  * Inquiries to: M. Houston, 36 Whetstone Clo. Farquhar Rd.         *
  4.  * Edgbaston, Birmingham B15 2QN ENGLAND                *
  5.  ********************************************************************
  6.  * File name: atan.c
  7.  * Program name: library modules only
  8.  * Source of file: The Public Domain Software Library.
  9.  * Purpose: maths function
  10.  * Changes: <who what when & why major changes have been made>      
  11.  ********************************************************************/
  12.  
  13.  
  14. /***********************************************************
  15.  *               The TULSA IBM C BOARD                     *
  16.  *                   918-664-8737                          *
  17.  *             300/1200 XMODEM, 24 Hours                   *
  18.  **********************************************************/
  19.  
  20.  
  21. #include "libc.h"
  22. #include "math.h"
  23. #include "errno.h"
  24.  
  25. static int nopper() {;}
  26.  
  27. #define PI              3.14159265358979323846
  28. #define PIov2   1.57079632679489661923
  29.  
  30. double atan2(v,u)
  31. double u,v;
  32. {
  33.         double f;
  34.         int (*save)();
  35.         extern int flterr;
  36.         extern int errno;
  37.  
  38.         if (u == 0.0) {
  39.                 if (v == 0.0) {
  40.                         errno = EDOM;
  41.                         return 0.0;
  42.                 }
  43.                 return PIov2;
  44.         }
  45.  
  46.         save = Sysvec[FLT_FAULT];
  47.         Sysvec[FLT_FAULT] = nopper;
  48.         flterr = 0;
  49.         f = v/u;
  50.         Sysvec[FLT_FAULT] = save;
  51.         if (flterr == 2)        /* overflow */
  52.                 f = PIov2;
  53.         else {
  54.                 if (flterr == 1)        /* underflow */
  55.                         f = 0.0;
  56.                 else
  57.                         f = atan(fabs(f));
  58.                 if (u < 0.0)
  59.                         f = PI - f;
  60.         }
  61.         if (v < 0.0)
  62.                 f = -f;
  63.         return f;
  64. }
  65.  
  66. #define P0 -0.13688768894191926929e+2
  67. #define P1 -0.20505855195861651981e+2
  68. #define P2 -0.84946240351320683534e+1
  69. #define P3 -0.83758299368150059274e+0
  70. #define Q0 +0.41066306682575781263e+2
  71. #define Q1 +0.86157349597130242515e+2
  72. #define Q2 +0.59578436142597344465e+2
  73. #define Q3 +0.15024001160028576121e+2
  74.  
  75. #define P(g) (((P3*g P2)*g P1)*g P0)
  76. #define Q(g) ((((g Q3)*g Q2)*g Q1)*g Q0)
  77.  
  78. double atan(x)
  79. double x;
  80. {
  81.         double f, r, g;
  82.         int n;
  83.         static double Avals[4] = {
  84.                 0.0,
  85.                 0.52359877559829887308,
  86.                 1.57079632679489661923,
  87.                 1.04719755119659774615
  88.         };
  89.  
  90.         n = 0;
  91.         f = fabs(x);
  92.         if (f > 1.0) {
  93.                 f = 1.0/f;
  94.                 n = 2;
  95.         }
  96.         if (f > 0.26794919243112270647) {
  97.                 f = (((0.73205080756887729353*f - 0.5) - 0.5) + f) /
  98.                                 (1.73205080756887729353 + f);
  99.                 ++n;
  100.         }
  101.         if (fabs(f) < 2.3e-10)
  102.                 r = f;
  103.         else {
  104.                 g = f*f;
  105.                 r = f + f *
  106.                         ((P(g)*g)
  107.                         /Q(g));
  108.         }
  109.         if (n > 1)
  110.                 r = -r;
  111.         r += Avals[n];
  112.         if (x < 0.0)
  113.                 r = -r;
  114.         return r;
  115. }
  116.  
  117.